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Two dimensional space-filling bearings are dense packings of disks that can rotate without slip. 
We consider the entire first family of bearings for loops of size four and propose a hierarchical con¬ 
struction of their contact network. We provide analytic expressions for the clustering coefficient and 
degree distribution, revealing bipartite scale-free behavior with tunable degree exponent depending 
on the bearing parameters. We also analyze their average shortest path and percolation properties. 


PACS numbers: 89.75.Hc, 89.75.Da, 45.70.-n 

I. INTRODUCTION 

Bearings are mechanical dissipative systems of rotors 
that, when perturbed, relax towards a bearing state, 
where all touching rotors rotate without slip. When these 
bearings cover the entire space they are called space¬ 
filling bearings [1]. Moreover,if such packings are sheared 
between moving surfaces, they can be used as a model 
to explain the existence of regions where tectonic plates 
can creep on each other for long periods of time without 
triggering earthquake activity, known as seismic gaps [2] . 
Space-filling bearings have also been used as a heuristic 
model for scale-free velocity fields, where the superdiffu¬ 
sion of massive particles can take place [3] . 

Herrmann et al. [1] presented a numerical algorithm to 
construct configurations of two dimensional space-filling 
bearings of polydisperse disks for loops of size four on a 
stripe geometry. They showed that two families of bear¬ 
ings can be obtained, where each configuration is clas¬ 
sified by two integer indices m and n. The contact net¬ 
work of a bearing is obtained by mapping it into a graph, 
where nodes are the disks and links are established be¬ 
tween touching disks. In the bearing state, which has 
no slip, two disks rolling on each other must have op¬ 
posite sense of rotation. The contact networks are thus 
bipartite, with the type of node defined by its sense of 
rotation. The topological properties of the contact net¬ 
work are intimately related to the force chains [4] and the 
dynamical response of the bearing to perturbations [5] . 

Andrade et al. [6] have shown that the contact net¬ 
work of Apollonian packings is a scale-free, small world. 
Euclidean, space-filling, and matching graph. The in¬ 
teresting properties of this network, named Apollonian 
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network, have motivated a series of follow-ups to study 
their geometrical [7], magnetic [8-10], spectral [11], and 
dynamical properties [12, 13]. Even an extension to ran¬ 
dom networks has been proposed [14]. In contrast to 
bearings where loops are necessarily of an even number of 
disks, for the Apollonian network loops are of size three. 

Here, we consider the first family of space-filling bear¬ 
ings in the stripe geometry and analyze their contact net¬ 
works. Doye and Massen looked at these networks in the 
limit m = n and provided heuristic arguments to esti¬ 
mate their degree exponent [7]. We propose a hierarchical 
construction of such networks which allows us to analyze 
the entire range of indices m and n and provide analytic 
expressions for the degree distribution and clustering co¬ 
efficient. We also describe several other properties. The 
paper is organized as follows. In Section H we start with 
the special case m = n = 0. The general case is dis¬ 
cussed in Section HI. We finally draw some conclusions 
in Section IV. 


II. NETWORK FOR m = n = 0 
A. The network construction 

We begin with the specific case of the space-filling 
bearing of m = n = 0 (see Eig. 1). This bearing has 
translational symmetry with a unit-cell composed of two 
topologically identical loops of four disks, defined by the 
largest disks, where the top and bottom surfaces are 
treated as disks of infinite radius. Eor m = n the bear¬ 
ing has also C2 rotation symmetry around the center of 
the common edge of the two largest loops. Thus, it is 
sufficient to consider the hierarchical construction rule 
for the contact network of one loop in the unit cell. By 
construction, all loops consist of an even number of disks 
and the network is bipartite, with two types of nodes 
denoted as a- and 6-nodes. The construction rule is sum¬ 
marized in Eig. 2. One starts with a loop arrangement 
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FIG. 1. (Color online) Section of the m,n = 0 bearing and 
the hrst generation of its contact network. 


of four nodes (two a- and two 6-nodes), corresponding 
to the four disks of the largest loop. An a-node is only 
connected to 6-nodes. The first generation ^ = 1 is con¬ 
structed by adding an a-node to the center of the loop 
and connecting it to the two 6-nodes, splitting the loop 
into two loops. This new a-node corresponds to the cen¬ 
tral disk touching the two lateral ones in each loop of the 
unit cell shown in Fig. 1. Inside each loop one 6-node is 
included and connected to the two closest a-nodes. These 
two new nodes correspond to the other two disks verti¬ 
cally aligned with the previous one. At the end, the 
initial square is divided into four loops. The next gener¬ 
ations are obtained hierarchically by repeating the same 
procedure inside each loop. By construction, the contact 
network is planar and self-similar. 


and 


= = , (4) 

t=l 

respectively. The total number of nodes N is 

N{g) = Na{g)^N,{g)=4^ -1 . (5) 

At each generation, all a- and 6-nodes receive one new 
link for each adjacent loop. Since the number of such 
loops equals the degree, the latter doubles at each gen¬ 
eration. The new a-nodes have degree four, while the 
6-nodes have degree two. Hence, at generation g, the de¬ 
gree k{t) of a node, added at generation go that is part 
of the network for t = ^ — generations is. 


fca(<7-50) = 4 - 29- 90 ^ 

( 6 ) 

fc6(ff-ffo) = 2 - 29 -90 . 

( 7 ) 


At generation g the degree of a node is related to the 
generation go at which the node was added. This gener¬ 
ation is given by. 


9 - goa{k) 


\n{k/2) 
In 2 


and 


g - gob{k) 


ln(fc/2) 
In 2 


(8) 


( 9 ) 


B. Degree distribution 

We now provide an analytic expression for the degree 
distribution P{k), where k is the node degree (number of 
touching disks). Let us start with the number of nodes 
N{g) at generation g, and neglect the first four nodes. 
One starts with one loop of two a- and two 6-nodes at 
generation zero. At each generation, each loop is divided 
into four. Thus, the final number of loops is 4^. For 
each loop in generation ^ — 1, one a- and two 6-nodes are 
added to obtain the generation g, so that the number of 
a-nodes Na changes from generation ^ — 1 to ^ as, 

AN,{g) = 49-1 

and the number of 6-nodes as. 


The number of nodes of degree k at generation g equals 
the number of nodes added at generation go{k), given by 
Eqs. (1) and (2). Thus, the degree distribution Pa{k,g) 
is. 


^ ANa{goa{k)) 
Na{g) 


= 3 


45 

49-1 



In the same way, Pb{k,g) is, 


^ ANtigobjk)) 

Nb{g) 

_ 3 4® 

“ 4 49 - 1 ( 2 ) 


( 10 ) 


( 11 ) 


AiVb(5) = 2-49-1 


(2) The total degree distribution P{k, g) is then, 


The number of nodes at generation g is then. 


7V„(5) = £4*-i = 


49 - 1 


t=i 


P{k,g) 


AjVa(g0a(fc)) + AiVfe(g06(fc)) 

N{9) 

3 49 fk'Y'^ 

2) ■ 


3 


( 3 ) 


249-1 


( 12 ) 
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two 6-nodes) in a loop arrangement, b) A new a-node is added to the center of the loop and connected to the 5-nodes on the top 
and bottom of the loop, dividing it into two loops, c) Two new 5-nodes are also added, one inside of each loop and connected 
to the two closest a-nodes. The initial loop is now divided into four equal loops, d) The new generation is constructed in the 
same way inside each loop. 


In the limit ^ ^ oc, 


Pa{k)=3( 

s ■ 

(13) 

CO 1 

II 

{T ■ 

(14) 

Pit) = f 

(7^ ^ 

(15) 


Thus, the degree distribution scales as P(/c) oc with 
7 = 2, corresponding to a scale-free network. This ex¬ 
ponent is larger than the one obtained for the Apollo¬ 
nian network, where 7 = ^ 1.585 [ 6 ]. Note that the 

a/5 asymmetry disappears when one considers the two 
topological identical loops in the entire unit cell, for the 
a-nodes in the top loop correspond to the 5-nodes in the 
bottom one. 



N 


C. Shortest paths and clustering coefficient 

Spatial, self-similar networks are expected to exhibit 
some form of small-world behavior due to the confine¬ 
ment of connections [7], which is, for example, the case 
of the Apollonian network [15]. Numerically, this can 
be checked by analyzing the size dependence of the av¬ 
erage shortest path /, defined as the average minimum 
number of links necessary to form a connecting path be¬ 
tween pairs of nodes in the network. Figure 3 shows that 
/ = alnA/" -b 5, where a and 5 are constants, correspond¬ 
ing to a logarithmic scaling, as expected for small-world 
networks [16]. 

Small-world networks are typically highly clus¬ 
tered [16]. To quantify the degree of clustering one mea¬ 
sures the clustering coefficient (7, defined for each node as 
the fraction of pairs of neighbors that are directly con¬ 
nected, forming a triangle. In the case of bearings, all 
loops have a even number of nodes, the contact network is 
bipartite, and two neighbors of a node are never directly 


FIG. 3. Size dependence of the average shortest path /, for 
the contact network of bearings of m = n = 0, where N is 
the number of nodes. The shortest path scales logarithmically 
with the number of nodes, I = 0.61 In N + 0.65 as expected 
for small-world networks. 


connected. Lind et al [17] proposed a new clustering co¬ 
efficient for bipartite networks (^ 4 , defined as the fraction 
of pairs that are indirectly connected through one single 
node, forming a loop. Then, for a node of degree k, 


C4 


7 ^(Indir. conn, between neighb.) 
k{k-l)/2 


(16) 


In the following we will use this definition. First, we cal¬ 
culate C^a/bit)^ the clustering coefficient of an a/5-node 
that was added to the network t generations before. At 
each iteration, the degree of every node is doubled, by 
adding new neighbors. Each new neighbor is connected 
via a new node to other two neighbors (see Fig. 2). Thus, 
the number of indirectly connected pairs of neighbors in- 
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creases at each generation by twice the node degree. Ev¬ 
ery new a-node has degree four and from its six different 
pairs of neighbors, five are indirectly connected. Every 
new 6 -node has degree two and its pair of neighbors is 
always indirectly connected. By summing over genera¬ 
tions, one gets for a-nodes. 


Ciait) = 


5 + EU2-fca(0 

ka{t){ka{t) - l)/2 
2 


and for 6 -nodes, 


= ^2- - 

2 4 ■ 2* - 1 


^ ^ ^ ’ kb{i) 

k,{t){k,{t) - l)/2 

= 3 ■ 2 -* - ^ 


2 - 2 ^ 


(17) 


(18) 


where we employed Eqs. ( 6 ) and (7). Once again, note 
that the a/b asymmetry is only observed when we solely 
consider one loop. Eor the entire stripe, the top and 
bottom loops of the unit cell are equivalent, but the a- 
nodes of the top loop are 6 -nodes of the bottom one. 
Thus, when the entire unit cell is considered, the network 
is completely symmetric with respect to a and 6 . Both 
coefficients tend to zero as (^ 4 ( 6 ) ^ 2“^ for t ^ 00 . The 
clustering of the entire network can be evaluated from 
the average over all nodes. 


^ ^ ANa{t)C4ai9 -t) + ANi,{t)Cib{g - t) 

= h ■ 

(19) 

We evaluate this sum numerically as shown in Eig. 4, for 
different number of nodes in the network, and obtain that 
(^4 ~ 0.8625 in the thermodynamic limit. 


D. Bond percolation 

We now study bond percolation on a network corre¬ 
sponding to a unit-cell of the stripe, consisting of two 
initial squares sharing one link (see Eig. 1 ). We focus 
on the existence of a spanning cluster between the two 
nodes representing the top and bottom surfaces, respec¬ 
tively. To compute the percolation threshold, we per¬ 
formed Monte Carlo simulations for different values of 
bond occupation probability p and network size N. We 
estimate the threshold Pc as the value of p at which the 
probability of spanning is 1 / 2 . We employed the bisection 
method and considered values of p that differ by 0 . 001 . 
Eigure 5 shows the size dependence of the estimated value 
of Pc 5 where one clearly sees that Pc vanishes in the ther¬ 
modynamic limit. Asymptotically, the decay follows a 
power-law Pc ^ , with z/ ^ 7. The same threshold 

is observed for the Apollonian network and other scale- 
free networks with 7 < 3. However, the convergence to 
the thermodynamic value is much slower here than for 
the Apollonian network, where u ^ 3 [ 6 ] . 
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FIG. 4. Clustering coefficient C 4 for the contact network of 
bearings of m = n = 0 as a function of the number of nodes 

N. 
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FIG. 5. Dependence o the percolation threshold pc on the 
network size N. The black line corresponds to the least- 
squares fit to the data of a power-law, pc ~ ]\^-fraci 2 iy^ with 
exponent z/ = 7 ± 1. 


III. GENERAL CASE 
A. The network construction 

We now consider the general case of the contact net¬ 
work for a bearing in the first family for loops of size four 
with any m and n. Figure 6 shows examples of bearings 
generated with different m and n, with the respective 
contact network on top. For all cases, the unit cell con¬ 
sists of loops of size four with the largest disks, including 
the top and bottom surfaces, respectively. However, the 
number of such loops varies with m and n and the rota¬ 
tion symmetry is broken for m ^ n. At each iteration the 
number of vertical and horizontal new loops constructed 
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FIG. 6. (Color online) Examples of space-filling bearings of the first family, for different values of m and n. 


inside each loop also depends on m and n, respectively. 
As summarized in Fig. 7, we first discuss how to deter¬ 
mine the initial number of loops in the unit cell (left 
panel) and proceed discussing how to hierarchically fill 
each loop (right panel). 

For all values of m and n, the unit cell of the bearing 
consists of a top and a bottom part. The number of 
initial loops on top and bottom equals n + 1 and m + 1, 
respectively (see loops of blue-square nodes in Fig. 6 and 
left panel in Fig. 7). Note that the (m,n) configuration 
is equivalent to the (n, m) configuration after a rotation 
of TT around the point where the common edge of the top 
and bottom loops crosses the middle of the stripe. To 
hierarchically construct the network one starts with the 
n + 1 top and m + 1 bottom loops. Hereafter, we solely 
consider one top and a bottom loop (sharing one edge), as 
the construction of the other loops is straightforward. To 
form the first generation we first add m + 1 a-nodes to the 
top loop and connect them to the two existing (lateral) 
6-nodes, dividing the initial loop into m + 2 loops. Then, 
in each new loop, we add n + 1 6-nodes and connect 
them to the top and bottom a-nodes. We are left with 
(n + 2)(m + 2) loops inside the top loop. Second, we 
construct the interior of the bottom loop. There, we 
start by adding n + 1 6-nodes and connect them to the 
two existing (lateral) a-nodes. Then, in each one of the 


new n+2 loops we add m+1 a-nodes and connect them to 
the top and bottom 6-nodes. In the bottom, we are also 
left with (n + 2)(m + 2) loops. Proceeding iteratively 
in the same way, we hierarchically construct the entire 
contact network of the space-filling bearing, for any m 
and n. 


B. Degree distribution 

We now provide an analytic expression for the degree 
distribution Pm,n{k) for any m and n, following the same 
strategy as for m = n = 0 in Sec. IIB. For simplicity, we 
restrict the calculation to networks inside a single loop 
(one initial top loop), as the total degree distribution of 
the unit cell is straightforwardly obtained as a weighted 
average of Pm,n{k) and Pn,m{k), corresponding to the 
degree distribution in the top and bottom loops, where 
the statistical weights are given by the initial fraction of 
top (n + 1) and bottom (m + 1) loops, respectively. 

At each generation, (n+2)(m+2) loops are constructed 
for every loop in the previous generation. As before, we 
neglect the initial set of nodes. At generation g — 1 there 
are (n + + 2)^“^ loops. Thus, from generation 

^ — 1 to ^f, the change in the number of a-nodes is, 

ANa{g) = {m+l)[{n + 2){m + 2)Y~^ , 


( 20 ) 
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n + 1 loops 



m + 1 loops 



FIG. 7. (Color online) The left panel shows a sketch of the initial network set up (blue, sites drawn as squares) and the first 
generation in the general case (red, sites as dots). In the right panel are the first two generations and part of the second (green) 
for m = 1 and n = 2. 


and in the number of 6-nodes is, 

ANb{g) = {m + 2){n + l)[{n + 2){m + 2)]’^~^ , (21) 

corresponding to {m + 1) new a- and {m + 2)(n + 1) new 
6-nodes per loop. The number of nodes at generation g 
is then, 


9 

Na{g) = + 1) 

t=l 

^ (m+l)([(n + 2)(TO + 2)]^-l) 

(n + 2)(m + 2)-l ’ ^ ’ 

and 


A 6-node has initially degree two and its degree increases 
by a factor of m + 2 at each generation. Thus, at gen¬ 
eration g, the degree of a 6-node added at generation go 
is, 


h{g-go) = 2{m + 2y-^^ . (26) 


Consequently, the node of degree k at generation g that 
was added at generation go is given by, 


g - goa{k) 


ln(fc/2) 
ln(n + 2) 


(27) 


for a-nodes and 


9 

Nb{g) = (m + 2)(n + 1) ^ [(^ + 2){m + 2)]^“^ 

t=i 

^ + 2)(n + l)([(n + 2)(m + 2)]^ - 1) 

(n + 2)(m + 2)-l ^ 

respectively. The total number of nodes N is, 

N{g) = Na{g) + Nb{g) = [{n + 2)(m + 2)]" - 1. (24) 

The degree /c of a node increases monotonically with 
the generation. An a-node has initially degree 2(n + 2) 
and its degree increases by a factor of n + 2 at each gen¬ 
eration. Thus, at generation g, the degree of an a-node 
added at generation go is, 


g - gob{k) 


ln(fc/2) 
ln(m + 2) ’ 


(28) 


for 6-nodes. The degree distribution for the a-nodes in 
the loop is then, 


Pm,nik,g) 


ANa{goa{k)) 

Naig) 

[(n + 2)(m + 2) — 1]/m,n(5) ( 2 ) 


where, 


[(n + 2)(m + 2)]® 
[(n + 2)(m + 2)]® - 1 


fca( 5 -fi'o) = 2(n + 2)® ®°+^ 


(25) 


fm ,n ig) 


(30) 
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For the 6-nodes is, 


iik,g) 


AN,{go,{k)) 

Nb{g) 

(n + 2 )(m + 2 ) — 1 


fm,n{9) ( ^ 


-[l+M=±^ 


[n + 2 ) {jn + 2 ) 

And the total degree distribution Pm,n{k^g) is, 

ANa{goa{k)) + AN,{gob{k)) 


ln(m + 2) J 

(31) 


iik,9) = 


— fm,ni9) 


N{9) 

( to + 1 ) 


\l I ln(m + 2) 1 

ln(n + 2) J 


Note that the degree exponent is symmetric to permu¬ 
tations of m and n. The degree distribution of the entire 
unit cell is symmetric at every generation as it is given 
by, 


p{k) = 


{n + l)Pm^n{k) + (m + l)Pn^m{k) 


m + n + 2 

[m-\- 
m-\-2 


\2l r 1 I ln(m + 2) ] 

(n+l)(m +1) + (1)”^ 


Tfl T 77/ T 2 


+ - 


\'2 r 1 _L ln(n + 2) I 

n+2 +(^ + l)(^+l) (I) i-(^+2)J 


777/ + 77/ + 2 


(35) 


r-. ln(rr + 2) 1 

77/ + 1 //C\ ^ ln(m + 2)J 

77/ T 2 \ 2 y 


(32) 


C. Shortest path and clustering coefficient 


In the limit ^ ^ oc, 


Pm,n{k) = {m+1) 


N ln(m + 2)l 


N ln(rr + 2) 1 

n + 1 

77/ T 2 \ 2 y 


(33) 


When 777 / = 77 /, 7 = 2 is recovered [7]. For m ^ the 
degree distribution is a sum of two power laws. Asymp¬ 
totically, for /c ^ oo, Pm,n{k) is dominated by the term 
with the smallest exponent and thus 


7 


= min 



\n{m + 2 ) 
ln(77/ + 2) 


ln(n + 2 ) ) 
\n{m + 2 ) J 


(34) 


The contact network of a bearing is always a scale-free 
network of 1 < 7 < 2, as shown in Fig. 8 . 


We numerically analyze the size dependence of the 
average shortest path I for all combinations of n^m = 
0,1, 2, 3,4. For all cases, we find a logarithmic scaling of 
I with the number of nodes, consistent with small-world 
networks. For m = n we find that the prefactor of the 
logarithmic scaling is independent on the value of the in¬ 
dices, as also observed for 7 . If m ^ then the prefactor 
changes with m and n as shown in Fig. 9, for example, 
for fixed ttT/ = 0 , the prefactor decreases with n. 

Next, we calculate (^ 4 . A new a-node has 
(n+3)(n+2)-l pairs of neighboring 6 -nodes connected in¬ 
directly through one a-node. At each generation, when 
new neighbors are added to each loop adjacent to this 
a-node, the number of connected pairs increases by 
(n+3)(n+2)/2-l. Hence, an a-node that was added to 
the network at generation ^ 0 , and is part of the network 
for t = ^ — ^0 generations, has clustering coefficient. 


C^a{t) = 


(tT/ + 3) (tT/ + 2) — 1 + 


(n-|-3)(n-|-2) —1 
2 


E*l7a(*) 


ka{t){ka{t) - l )/2 

(n + 3 )(n + 2 )-l 2 (n + 2 ) 


{n + 2)^+1 [ 2 ( 77 / + 2)^+1 - 1] 
(tT/ T 2) ^ as t —y cxd . 


(n+3)(n+2)-l 

2 


- 1 


(n+2)^-l 

n+l 


{n + 2)^+1 [ 2 ( 77 / + 2)^+1 - 1] 


Initially, for 6 -nodes there is only one pair of neighbors added per adjacent loop at each iteration. Thus, 
indirectly connected and (^+^K^+^) _ i connections are 


C^b{t) = 


1 + 


(m-|-3)(m-|-2) —1 
2 


- 1 




kb{t){kb{t) - l )/2 


1 + 2 


(7ti-|-3) (?ti-|-2) — 1 


- 1 


(?ti-|- 2) ^ — 1 
m+l 


(m + 2)*(2(m + 2)* - 1) 
{m + 2 )“* as t —> oo. 


r\j 



















































FIG. 8 . (Color online) Left: Degree exponent 7 as a function of m and n. Note that, 7 = 2 for m — n. Right: Degree 
exponent 7 as a function of n for m = 100 . 
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FIG. 9. Left: Average shortest path I as a function of the number of nodes N for n = m. Lines were shifted vertically by 
adding n for better visualization. Right: Average shortest path I as a function of the number of sites N for m = 0. Lines were 
also shifted vertically. The lines are just guides to the eye. 


The argument of the power law is different for a- and b- D. Bond percolation 

nodes as it depends on n and m, respectively. Both C 4 a 
and asymptotically vanish. The faster the degree of 

a node type grows the faster its C 4 falls off. We performed simulations of the bond percolation 

model on a unit-cell for different pairs of indices m, n. As 
in Sec. IID, we define the spanning cluster as a set of con¬ 
nected nodes that include the top and bottom surfaces. 
Note that the top and bottom surfaces correspond to dif¬ 
ferent types of nodes, a and 6, respectively (see Fig. 6). 
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For all considered values of m and n we find that 
the percolation threshold Pc vanishes in the thermody¬ 
namic limit (infinite system size) and the estimator for 
the threshold scales as Pc{N) ^ N~^. Our results for 
m = n suggest that u does not change with the bearing 
indices {m and n), like we also found for the degree ex¬ 
ponent 7 in Sec. IIIB. Since the number of nodes in the 
network grows exponentially with the generation, we re¬ 
frain from performing a detailed size-dependence analysis 
to obtain u with high precision. 

IV. FINAL REMARKS 

We studied the contact network of space-filling bear¬ 
ings of loops of size four in the first family. We pro¬ 
posed a hierarchical rule to construct the network and 
provided analytic expressions for the degree distribution 
and clustering coefficient. We also studied numerically 
the shortest path and percolation properties. We showed 
that the exponent 7 changes in the range ( 1 , 2 ] and that 
is always two when m = n. Numerical simulations also 
suggest that the correlation exponent u for the percola¬ 
tion transition does not change with the bearing indices. 
Our networks are bipartite and we find that if m ^ n the 
degree distribution of the two species scale with differ¬ 


ent exponents inside each loop. To our knowledge, this 
is the first example of an artificial hierarchical network 
exhibiting this property which has already been observed 
empirically for sexual networks [18]. 

We are proposing a method to generate deterministic 
hierarchical scale-free networks of different 7 exponents, 
which are amenable to analytic treatment. As it was 
accomplished for the Apollonian network, possible ex¬ 
tensions of our work include the study of their magnetic, 
spectral, and dynamical properties [19-21]. Other pos¬ 
sibilities include the study of the networks of random 
space-filling or three-dimensional bearings, larger loops 
and second family of bearings. 
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